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Driving rate dependence of avalanche statistics and shapes at the yielding transition 
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We study stress time series caused by plastic avalanches in athermally sheared disordered mate¬ 
rials. Using particle-based simulations and a mesoscopic elasto-plastic model, we analyze size and 
shear-rate dependence of the stress-drop durations and size distributions together with their aver¬ 
age temporal shape. We find critical exponents different from mean-field predictions, and a clear 
asymmetry for individual avalanches. We probe scaling relations for the rate dependency of the 
dynamics and we report a crossover towards mean-field results for strong driving. 
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Many materials respond to slow driving with strongly 
intermittent dynamics. Examples include Barkhausen 
noise in ferromagnets [1-3], stick-slip motion in earth¬ 
quakes [4], serration dynamics in plasticity of solids [5], 
and avalanche dynamics in crack propagation [6, 7], 
driven foams [8] and domain wall motion [9]. 

As in equilibrium critical phenomena, global quantities 
linked to such bursting collective events are usually power 
law distributed and allow for the introduction of scaling 
functions. In the slow driving limit, the onset of motion 
can be interpreted as an out-of-equilibriunr phase transi¬ 
tion, suggesting the existence of families of systems that 
display similar avalanche statistics. To better identify 
this universality classes, both experimental [10-17] and 
theoretical [13, 18-21] works have discussed the avalanche 
“shapes”, going beyond the study of scaling exponents. 

In deformation experiments of amorphous systems, 
such as grains, foams or metallic glasses, avalanche dy¬ 
namics are typically evidenced in the time series of the 
deviatoric component of the stress tensor. In the limit 
of vanishing deformation rate we approach the so-called 
“yielding transition”. The question whether yielding can 
be characterized as a continuous dynamical phase tran¬ 
sition, belonging to a specific universality class, is still 
under debate. The analysis of avalanche statistics close 
to yielding has therefore a particular relevance. 

In this letter, we study the emerging yielding dy¬ 
namics in a simple shear geometry with imposed driv¬ 
ing rate. Our focus lies on the shear-rate dependence 
of the avalanche statistics and thus complement recent 
quasi-static studies [22-25] , To address the low shear-rate 
regime we use a coarse-graining approach, proven to yield 
qualitative and quantitative relevant predictions [26-31], 
and compare the low shear-rate results of our meso-scale 
model with quasistatic particle-based simulations. 
Molecular dynamics (MD) - We consider a mixture of A 
and B particles interacting via a Lennard-Jones poten¬ 
tial: 14 b (r) = 4e AB [(a AB /r) 12 - (cr AB /r) 6 ] with r being the 
distance between two particles. Units of energy, length 
and mass are defined by e AA , cr AA and m A ; the unit of 
time is given by r 0 = cr AA \J(m A /e AA ). The potential is 


truncated at R c = 2.5 and a force smoothing is applied 
between an inner cut-off Ri n = 2.2 and R c . The two 
species of particles have equal mass to, but different in¬ 
teraction parameters to prevent crystallization. We set 
£aa = 1-0, £ab = 1-5, e BB = 0.5, cr AA = 1.0, a AB = 0.8 
and (t B b = 0.88 and m = 1. The ratio of particles of 
species A and B is chosen N A /N B = 13/7 and 8/2 for 2d 
and 3 d systems, respectively. Glassy states are obtained 
(with LAMMPS [32]) by quenching to zero temperature 
at constant volume systems equilibrated at T = 1. An 
athermal system is achieved by applying to each particle 
a viscous drag force F drag = —Tv, where v is the par¬ 
ticle peculiar velocity. We condition the dynamics to be 
strongly overdamped [22, 33] (T = 1). Avalanche statis¬ 
tics are obtained following a quasistatic protocol [22, 23]. 
We impose simple shear at rate 7 = 10~ 6 by deforming 
the box dimensions and remapping the particle positions. 
Following [22] , the shear-rate 7 is set to zero when a steep 
increase in kinetic energy occurs (onset of plastic defor¬ 
mation) and only restored when the kinetic energy drops 
below a threshold. 

Elasto-plastic (EP) model - We coarse-grain an amor¬ 
phous medium onto a mesoscopic lattice: each node rep¬ 
resents a block of material holding exactly one shear 
transformation [33-36], for which we assume the same 
geometry as the globally applied simple shear. To each 
site i we associate a local scalar shear stress a t and a 
state variable n*, indicating whether the site plastically 
deforms (n = 1) or not (n = 0). Local stresses evolve 
with the overdamped dynamics: 

dtcri = ni + ^ X (1) 

3 

with p, = 1 the elastic modulus, 7 the externally ap¬ 
plied shear-rate, t = 1 a mechanical relaxation time and 
dt'if = the strain rate produced by a plastic re¬ 

arrangement at site j. Gij denotes the discretized Es- 
helby propagator [37] , that obeys a quadrupolar symme¬ 
try in the shear plane with a dipolar long-range character, 
G(r,r') = cos(49 r r')/\r—r'\ d . A site yields (m = 0 —» 1) 
when its stress reaches a local threshold cr, > erf, and re- 
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FIG. 1. Stress-drop size distributions. Main panels show 
rescaled distributions L d f T Ps vs. S/L d i of the EP model 
compared to MD quasistatic simulations (arbitrary shift ap¬ 
plied for the comparison). Insets show not-scaled curves, (a) 
3dEP model data for linear system-sizes L = 16 (green cir¬ 
cles), 32 (orange squares), 64 (blue diamonds), 128 (plum 
triangles) and shear-rate 10 ~ 4 (full symbols). For L = 32, 
7 = 10 4 , 10~ 5 are also shown (light and dark orange open 
squares). Gray scale triangles correspond to quasistatic 3dMD 
with L — 40,60,80 (from light to dark), (b) 2dEP data 
for linear system-sizes L = 256 (green circles), 512 (orange 
squares), 1024 (blue diamonds) and 2048 (plum triangles) 
at 7 = 10 -5 . Gray scale triangles correspond to quasistatic 
2dMD with L = 80,160, 320 (from light to dark). 


covers its elastic state (rii = 1 —t 0 ) when a prescribed 
local deformation increment is attained after yielding, 
f\d t (Ji/n + d t ”ff\dt > 7 C . Each time a site yields a 
new yield stress erf is drawn from a distribution of mean 
a o- Model details and parameter choices can be found in 
Ref. [38] and in the Supplemental Material [39]. 
Stress-drop statistics and shear-rate dependence- From 
the stress-time series we individualize stress-drops, and 
define an extensive quantity S proportional to the abso¬ 
lute stress difference multiplied by the system volume. 
We compare in Fig.l the stress-drop distributions Ps 
in the limit of low 7 for the EP model with the qua¬ 
sistatic MD results. In both two (2d) and three dimen¬ 
sions (3d), apart from a plateau regime for small stress- 
drops that depends on shear-rate, numerical integration 
step and system size, we fit the data using a power-law 
P s ~ S~ T f(S/S c ), with / an exponentially decaying cut¬ 
off function [40] (exponent definitions in Table I). Notic¬ 
ing that the distributions Ps become independent of 7 
in the zero shear-rate limit and in agreement with previ¬ 
ous works [23, 25], we postulate a system size dependent 
cut-off S c ~ L dt , with df the fractal dimension of the 
avalanches [23, 25, 41]. The comparison of these stress- 
drop statistics with MD results reveals a fair agreement, 
up to an arbitrary scaling factor related to the difference 
in simulated length scales. 

The fitted values of r for the EP model, both in two 
and three dimensions (r 2d ~ 1.28, r 3d ~ 1-25), com¬ 
pare very well with our and earlier obtained MD re¬ 
sults [22, 23], are compatible with previous lattice mod- 


FIG. 2. Shear-rate dependency of the dynamics for the SdEP 
model, (a) Log-log plot of Aoo = (a — o c )/oo vs 7 . Circles 
correspond to the best estimation of (t c /oq = 0.687; along¬ 
side lines, to choices of 0.683 and 0.691 instead. Full and 
dashed lines are power-law fits in selected ranges (extrapo¬ 
lated for comparison). Inset: Crossover of 1 /0 as explained 
in the text, (b) Steady-state distributions P x of the local 
distances to threshold x = o y — <7 for different shear-rates 
7 £ {10 -1 ' 4 ,..., 10“°}. Inset: Stress-drop distributions for 
7 £ { 10 “ 1 ,..., 10 ~ 3 }, rescaled and shifted as explained in 
the text. Arrows indicate the sense of increasing shear-rate, 
(c) Rescaled distributions of stress-drop duration 7 _aT P T vs. 
T'y 01 for 7 = 10 - 2 ,10“ 3 , 10 - 4 ,10 ~ 5 (from light blue to dark 
plum, left to right in inset), and system-sizes L = 64 (closed 
symbols) and 128 (open symbols). The dashed line shows a 
law Pt ~ T -1 ' 44 . X nS et: Unsealed data, (d) Average size 
S for stress-drops of the same duration as a function of T 7 “ 
for L = 64 and 7 = 10“ 2 ,10~ 3 ,10~ 4 ,10 -5 . The dashed line 
shows S ~ T 1,58 . Inset: Unsealed data, shear-rate decreases 
from left to right. 


els [42], and lie within error bars of those provided by 
FEM models [43]. Still, they disagree with what was 
obtained with quasistatic protocols in cellular automa¬ 
ton models [25] (especially in 3 d where r 3d ~ 1-43), and 
they contrast even more with the usual mean-field pre¬ 
diction [44] r MF = 3/2 (see [45] for an alternative analy¬ 
sis). The values obtained for df (d 2d ~ 0.9, e? 3d ~ 1.3) 
are compatible with quasistatic MD simulations, but 
slightly smaller than those reported in automaton mod¬ 
els [25]. They suggest a line geometry of the correlated 
slip events [24, 46], with a modest but clear trend towards 
a more compact structure in 3d. 

Some main results concerning the finite driving rate are 
summarized in Fig. 2 for the 3<£EP model, similar results 
are found for the 2d case (not shown). The consequences 
of applying a finite shear-rate are twofold[47]: 

(I) The first important observation is that with increas- 
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ing driving rate the critical exponents tend towards the 
mean-held predictions. The yielding exponent /? for ex¬ 
ample, defined through 7 oc (cr — a c )P , can be derived 
from the fits in Fig.2(a) rendering a non-trivial value 
/3 ~ 1.55 in the low shear-rate regime. For larger shear- 
rates this value crosses over to /3 ~ 2 predicted by the 
Hebraud-Lequeux model [48] . By sliding a fixed size log¬ 
arithmic window in 7 (comprising ^12 points of the main 
plot data set) and fitting within, we show the resulting 
1//3 as a function of the starting position of the window 
in the inset of Fig. 2(a). Similarly we observe a crossover 
of the exponents in the steady-state distribution P x of 
the local stress excess [49, 50] 27 = erf — 07 , Fig.2(b). 
Again in the limit of vanishing shear-rates we observe 
the curves approaching a shape that initially grows as 
P x ~ x s with a non-trivial exponent, as found in the qua¬ 
sistatic case [25, 49], attributed to an anomalous random 
walk process of the local stress with an absorbing bound¬ 
ary condition at x = 0 [51]. Flowever, as we increase 
the shear-rate P x changes, eventually yielding 6 ~ 0. 
The driving progressively dominates over the signed kicks 
from elastic interactions, yielding a biased diffusion of the 
x’s values. This ultimately produces a strictly positive 
local stress evolution, resembling the x dynamics of the 
depinning problem [25]. The inset of Fig.2(b) shows a 
feature compatible with the shear-rate dependence of P x 
and with the /3 crossover. For different shear-rates, we 
plot k-yS 15 Ps vs. S, where is an arbitrary scaling co¬ 
efficient to separate the curves and improve visualization. 
We observe a range of low shear-rates where the slope of 
the transformed distributions is almost unchanged and 
fully consistent with Fig. 1(a). Above a rate of deforma¬ 
tion of about ~ 0.015, curves progressively flatten, even¬ 
tually becoming horizontal. Plotting S 15 Ps , we show the 
departure of Ps from the MF expectation Pg F oc S -15 
as the critical point is approached. When investigating 
the distribution of stress fluctuations 77 = Gij ^ 

on each site, we find consistently a change from a peaked 
distribution with fat tails towards Gaussian-like distribu¬ 
tions as we increase the shear-rate. We infer from this, 
that the strong correlations at vanishing shear-rates (rea¬ 
son for the non-trivial criticality) become negligible for 
stronger driving, so that the exponents end up being well 
described by mean-field assumptions. 

(II) The second consequence of a finite driving rate is 
that the critical scaling regime shows not only finite size, 
but also finite shear-rate effects [46, 52]. When imposing 
a finite deformation rate, each stress-drop is character¬ 
ized not only by its magnitude or size 5, but also by its 
duration T. For each stress-drop we define a given dura¬ 
tion T, as the time elapsed between the beginning and the 
end of the drop. In Fig.2(c) we present the distributions 
of durations Pt for a fixed system size and different shear- 
rates. In the probed shear-rate regime we find the depen¬ 
dence on L to be negligible, thus Pt(T , L , 7 ) = Pp(T, 7 ). 
The main panel shows rescaled curves assuming the func- 



FIG. 3. Stress-drop shape properties for a 3dEP model-, (a) 
Rescaled stress-drop shape Vr(t) = Vr(t)/ maxt(Vr(t)) av¬ 
eraged over stress-drops of duration Tie, as a function of 
rescaled time t = t/T. From left to right, we show curves at 
increasing T. Inset: Bordered squares represent fitted values 
of the asymmetry parameter a g for different choices of ( 7 , L). 
Color code depicts the fit a g = 10 - 0 ' 42 T - 0 ' 43 7 _ 0 ' 37 L -1 ' 25 for 
T = 0.5. (b) Amplitude B of the stress-drops vs. T, for 
L = 32 and 7 = 10 -4 , as obtained from the fits. The inset 
shows corresponding a s (T) and a g (T) (see text). 

tional dependence Pt ~ T~ T glT-y 01 ), with g an expo¬ 
nentially decaying function. We obtain for the 3 d case, 
T 3d = 1.44 and a 3d = 0.3. Naturally, we expect the scal¬ 
ing of Pp to be dominated by a growing length scale £ in 
the critical limit, where the relations T ~ £ z and S ~ £ df 
hold. Therefore, we expect a scaling relation S ~ T s 
with 6 = df/z, that we observe over a range of shear- 
rates, yielding the exponent S 3d ~ 1.58 (see Fig.2(d)), in 
contrast with the mean-held <5 M f = 2. More generally, 
we observe empirically a power-law scaling of S with T, 
7 and L. Actually, extending the dependencies of the 
cut-off values in size, L df , and duration, 7 ““, the mean 
S at each T should follow S(T,L, 7 ) = C(L,"/)T S with 
C(L, 7 ) ~ L df 7 a<5 . This relation is fairly verified for the 
dependence on 7 , illustrated in Fig.2(d). A rescaling of 
the size dependence leads to an exponent larger by 15% 
than df estimated from Pg. 

Stress-drop shapes- We address now the analysis of the 
functional form of the stress-drops, i.e., the time evolu¬ 
tion of the stress-drop velocity [15-17, 44]. In Fig.3(a) 
we show rescaled stress-drop velocities Vr (stress-drop 
shapes) for a 3d system, averaged over drops of the 
same duration T within the power-law scaling regime of 
Fig.2(d). We observe that drops of short duration show 
a noticable asymmetric shape, with faster velocities at 
earlier times. As duration increases, the shape becomes 
gradually more symmetric. To analyze this asymmetry 
of stress-drop shapes for different durations, system sizes 
and applied shear-rates, we fit them with a formula pro¬ 
posed in Ref.[16] Vr{t) oc — f)) c (l — a s (t — 0.5)) (see 
also [15, 53]), with B the amplitude of the shape and 
a s a parameter quantifying the deviation from a sym¬ 
metric inverted parabola. We confirm the expected rela¬ 
tion c = S — 1 (recall S oc T s and compare Fig.3(b)and 
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Fig.2(d)). In our range of parameters c is almost in¬ 
dependent of L and 7 . More relevant for our analysis 
is the behavior of the fitting parameter a s (see Fig.3(b) 
inset), that shows clearly the crossover from nearly sym¬ 
metric to asymmetric shapes as we focus on shorter du¬ 
rations T. To avoid a fit with various parameters, we use 
an alternative, purely geometrical measurement of the 
asymmetry that is relevant even beyond scaling regime, 
a 9 = fo di - When com Puting a g (T) for 

different shear-rates at fixed T and L, a g increases as 7 
decreases; whereas for fixed T and 7 , a g decreases as L 
increases (see inset of Fig.3(a)). In the quasistatic limit, 
where just one independent avalanche occurs at a time 
we expect asymmetric stress-drop shapes characterizing 
individual avalanches. When we increase the driving rate 
at fixed system size or, equivalently increase the system 
size at a fixed rate, we expect stress-drops to result from 
many independent avalanches, since the density of plas¬ 
tic regions is determined and increased by the driving 
strength [46]. Here, the resulting stress-drop shape draws 
closer to the mean-field symmetric shape. 

Conclusions- We studied with a mesoscopic model the 
avalanche statistics close to the yielding transition, ver¬ 
ifying the relevance of our approach by comparing with 
particle-based quasistatic simulations. In Table I we 
summarize the critical exponents obtained for 2d and 3 d. 
Our results clearly reinforce the idea of a non-trivial uni¬ 
versality class for the yielding transition, in agreement 
with earlier findings [23, 25, 42]. Our estimated expo¬ 
nents, confirm within error bars the scaling relations pro¬ 
posed by Lin et al. [25] . We also note that our values of 
r and t' are indistinguishable from the exponents ex¬ 
pected for the Id long-range (1 /r 2 ) depinning universal¬ 
ity class [ 6 , 7]. Although the loading path dependence of 
the critical exponents remains an open issue, this is an 
interesting accordance and points towards the role played 
by the avalanche slip-line geometry. 

In the regime of larger shear-rates we find that sev¬ 
eral exponents of the stress-drop statistics draw closer to 
mean-held predictions. The rise of an increasing num¬ 
ber of independent regions with yielding activity (paral¬ 
lel occurring avalanches) justifies the crossover to triv¬ 
ially random statistics. In particular our data reveals 
a yielding exponent approaching the prediction of the 
Hebraud-Lequeux model [48, 54, 55]. Further the finite 
shear-rate protocol allows for the introduction of an ad¬ 
ditional exponent a that should enter the scaling rela¬ 
tions. If we assume a usual scaling scenario, we expect a 
diverging length scale depending on the distance to the 
yielding point £ ~ {a — a c )~ v , such that £ ~ ) since 

7 ~ (<r — a c )P. Then T ~ £ z yields directly the scaling 
relation a = zv/[3. We have not measured v , but assum¬ 
ing v = 1 /(d — df) [25] to be valid we get ct 2 d = 0.34 and 
ot?,d = 0.31, close to the measured values. 

Within the scaling regime for T we observe both asym- 



Expression 

This work (2d | 3d) 

^r-depinning Id 

MF 

p 

7 ~ (Act)^ 

1.54(2) 1.55(2) 

0.625(5) [56] 

2 [48] 

T 

Ps ~ S~ T 

1.28(5) 1.25(5) 

1.25(5) [6,7] 

1.5 [44] 

df 

S c ~ L d f 

0.90(7) 1.3(1) 

~ 1.38 [56] 

— 

/ 

T 

p T ~ T~ t ' 

1.41(4) 1.44(4) 

- 1.43 [ 6 ] 

2 [44] 

a 

T c — 7 “ 

0.38(4) 0.30(4) 

— 

— 

z 

T ~ C 

~ 0.57 ~ 0.82 

0.77(1) [56] 

— 

5 

S ~T S 

1.58(7) 1.58(5) 

~1.7 [ 6 ] 

2 [44] 

e 

P x ~x 6 

0.52(3) 0.37(5) 

0 

1 [49] 


TABLE I. Measured exponents for the avalanche statistics. 


metric and symmetric stress-drop shapes depending on 
system size, shear-rate and duration. This is why we 
propose to distinguish between individual avalanches (re¬ 
sulting from correlated plastic events) and stress-drop 
shapes (resulting from many independently occurring 
avalanches). 

The combined study of avalanche size and duration 
distributions and avalanche shapes has played an essen¬ 
tial role in our understanding of the universal aspects of 
crackling noise and depinning dynamics. With this work, 
we provide a first numerical prediction of similar quan¬ 
tities in the case of the yielding transition, with a clear 
indication of a complex non mean-field behavior. We 
hope this work will stimulate and provide a benchmark 
for future experimental studies on systems undergoing a 
continuous yielding transition, for which detailed data on 
noise statistics is presently very scarce. 
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